

%% data

dChar = readtable('../raw_data/district_characteristics.csv');
e09 = readtable('../raw_data/election_results_2009.csv');
e12 = readtable('../raw_data/election_results_2012.csv');
cs09 = readtable('../raw_data/campaign_spending_2009.csv');
cs12 = readtable('../raw_data/campaign_spending_2012.csv');

% order observations by candidate-coalition choice
d0 = find(e12.candCM==0); N0 = size(d0,1);
d1 = find(e12.candCM==1); N1 = size(d1,1);
d2 = find(e12.candCM==2); N2 = size(d2,1); N = N0+N1+N2;
dChar = dChar([d0;d1;d2],:);
e09 = e09([d0;d1;d2],:); e12 = e12([d0;d1;d2],:);
cs09 = cs09([d0;d1;d2],:); cs12 = cs12([d0;d1;d2],:);
clear d0 d1 d2

% district indicator/aggregation matrix
A = [1:N,1:N,1:N0,N0+N1+(1:N2),1:N0+N1,1:N]';
A = 1 * (A==A');

% lagged vote shares
P = [e09.MP, ...
    e09.NA, ...
    e09.PVEM + 0.5 * e09.PM, ...
    e09.PRI + 0.5 * e09.PM, ...
    e09.PAN] ./ e09.registered;

% 2012 vote shares
S = [e12.MP ./ e12.registered; ...
    e12.NA ./ e12.registered; ...
    e12.PVEM(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+N1+(1:N2),1) + e12.PVEM(N0+N1+(1:N2),1) + e12.CM(N0+N1+(1:N2),1)) ./ e12.registered(N0+N1+(1:N2),1); ...
    e12.PRI(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+(1:N1),1) + e12.PVEM(N0+(1:N1),1) + e12.CM(N0+(1:N1),1)) ./ e12.registered(N0+(1:N1),1); ...
    e12.PAN ./ e12.registered];


%% baseline specification: OLS

% demographics
D = [dChar.femaleHH,dChar.over60,dChar.rural];
dnames = {'female','over60','rural'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);
NN = size(X,1);

% add log-lagged vote shares
PP = zeros(NN,3);
PP([1:N0,N+(1:N0),2*N+(1:N0),2*N+N0+N2+(1:N0),3*N+N0+(1:N0)],1) = reshape(log(P(1:N0,:)),5*N0,1);
PP([N0+(1:N1),N+N0+(1:N1),2*N+N0+N2+N0+(1:N1),3*N+N0+N0+(1:N1)],2) = reshape(log([P(N0+(1:N1),1:2),sum(P(N0+(1:N1),3:4),2),P(N0+(1:N1),5)]),4*N1,1);
PP([N0+N1+(1:N2),N+N0+N1+(1:N2),2*N+N0+(1:N2),3*N+N0+N0+N1+(1:N2)],3) = reshape(log([P(N0+N1+(1:N2),1:2),sum(P(N0+N1+(1:N2),3:4),2),P(N0+N1+(1:N2),5)]),4*N2,1);
X = [X,sum(PP,2)];
K0 = size(X,2);
X0 = X; % save preferred specification for campaign stage

% (endogenous) campaign spending
C = [cs12.MP;...
    cs12.NA;...
    cs12.PVEM(1:N0,1);cs12.CM(N0+N1+(1:N2),1);...
    cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);...
    cs12.PAN];
XC = [X,C,C.^2];
K = size(XC,2);

% OLS
Y = log(S) - log(1-A*S);
Bols = XC \ Y;
xi = Y - XC * Bols;
SS = XC' * diag(xi.*xi) * XC;
VBols = (XC' * XC) \ SS / (XC' * XC); % heteroskedasticity-robust variance
seBols = sqrt(diag(VBols)); % standard errors
rBols = [Bols,seBols,2*(1-normcdf(abs(Bols./seBols)))];

disp(' ')
disp('for Table 4 (column (I)):')
disp(print_results(round(rBols,3),dnames,[],0,[],0))


%% baseline specification: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];
X2S0 = X2S;

% OLS
Y2S = [log(e12.PVEM(N0+1:end)) - log(e12.CM(N0+1:end));...
    log(e12.PRI(N0+1:end)) - log(e12.CM(N0+1:end))];
B2S = X2S \ Y2S;
B2S0 = B2S;
xi2S = Y2S - X2S * B2S;
xi2S0 = xi2S;
VB2S = (X2S' * X2S) \ ...
    (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
VB2S0 = VB2S;
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table A4 (column (I)):')
disp(print_results(round(rB2S,3),dnames,[],0,[],1))


%% baseline specification: BLP

% lagged-spending instruments
Lagged = [cs09.MP;...
    cs09.NA;...
    cs09.PVEM(1:N0,1)+cs09.PM(1:N0,1);cs09.PRI(N0+N1+(1:N2),1)+cs09.PVEM(N0+N1+(1:N2),1)+cs09.PM(N0+N1+(1:N2),1);...
    cs09.PRI(1:N0,1)+cs09.PM(1:N0,1);cs09.PRI(N0+(1:N1),1)+cs09.PVEM(N0+(1:N1),1)+cs09.PM(N0+(1:N1),1);...
    cs09.PAN];
Z = [X,Lagged,Lagged.^2];

% Gandhi-Houde nonlinear instruments
Z = GH_Instr([X(:,end),C,C.^2],Z,A);
Z0 = Z;
KZ = size(Z,2);

% sparse-grid integration nodes and weights (from Heiss and Winschel, 2008)
[nu,wnu] = nwspgr('KPN',2,6);

% first iteration
W = Z' * diag(xi.*xi) * Z; % GMM weighting matrix
Bblp = [Bols;0.25;0.15]; % starting values
xi = xi_BLP(Bblp,xi,S,XC,A,nu,wnu); % compute unobserved valence
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)]; % Jacobian sparsity pattern
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ))); % Hessian sparsity pattern
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu)); % options for Knitro
[B,~,exit1] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt'); % to check derivatives against finite differences, use 'knitro0.opt'
xi = B(length(Bblp)+(1:length(S)),1);

% second iteration
W = Z' * diag(xi.*xi) * Z; % GMM weighting matrix
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit2] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[B(1:length(Bblp),1);xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);
xi0 = xi;

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = -[XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z' * Ghat;
Ghat0 = Ghat;
VBblp = inv(Ghat'*((Z'*diag(xi.*xi)*Z)\Ghat)); % robust variance
VB0 = VBblp;
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
B0 = Bblp;
rBblp = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table 4 (column (II)):')
disp(print_results(round(rBblp,3),dnames,[],0,1,0))


%% baseline specification: BLP + neighbors' spending

% neighbors' spending and instrument
nb = neighb([cs12.MP,...
    cs12.NA,...
    [cs12.PVEM(1:N0,1);zeros(N1,1);cs12.CM(N0+N1+(1:N2),1)],...
    [cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);zeros(N2,1)],...
    cs12.PAN]);
X = [X,[nb(:,1);nb(:,2);nb(1:N0,3);nb(N0+N1+(1:N2),3);nb(1:N0+N1,4);nb(:,5)]];
XC = [X,C,C.^2];
K = size(XC,2);
lnb = neighb([cs09.MP,...
    cs09.NA,...
    [cs09.PVEM(1:N0,1)+cs09.PM(1:N0,1);zeros(N1,1);cs09.PRI(N0+N1+(1:N2),1)+cs09.PVEM(N0+N1+(1:N2),1)+cs09.PM(N0+N1+(1:N2),1)],...
    [cs09.PRI(1:N0,1)+cs09.PM(1:N0,1);cs09.PRI(N0+(1:N1),1)+cs09.PVEM(N0+(1:N1),1)+cs09.PM(N0+(1:N1),1);zeros(N2,1)],...
    cs09.PAN]);
Z = [Z,[lnb(:,1);lnb(:,2);lnb(1:N0,3);lnb(N0+N1+(1:N2),3);lnb(1:N0+N1,4);lnb(:,5)]];
KZ = size(Z,2);

% BLP
W = Z' * diag(xi.*xi) * Z; % GMM weighting matrix
Bblp = [Bblp(1:K0);0;Bblp(K0+1:end)];
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_nb] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = -[XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z' * Ghat;
VBblp = inv(Ghat'*((Z'*diag(xi.*xi)*Z)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_nb = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table 4 (column (III)):')
disp(print_results(round(rBblp_nb,3),dnames,[],1,1,0))


%% region FE: OLS

% demographics
D = [dChar.region<=2,dChar.region==3,dChar.femaleHH,dChar.over60,dChar.rural];
rfe = {'north','southeast'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);

% add log-lagged vote shares
X = [X,sum(PP,2)];
K0 = size(X,2);

% (endogenous) campaign spending
XC = [X,C,C.^2];
K = size(XC,2);

% OLS
Bols = XC \ Y;
xi = Y - XC * Bols;
SS = XC' * diag(xi.*xi) * XC;
VBols = (XC' * XC) \ SS / (XC' * XC); % heteroskedasticity-robust variance
seBols = sqrt(diag(VBols)); % standard errors
rBols_RFE = [Bols,seBols,2*(1-normcdf(abs(Bols./seBols)))];

disp(' ')
disp('for Table 4 (column (IV)):')
disp(print_results(round(rBols_RFE,3),dnames,rfe,0,[],0))


%% region FE: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];

% OLS
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S_RFE = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table A4 (column (II)):')
disp(print_results(round(rB2S_RFE,3),dnames,rfe,0,[],1))


%% region FE: BLP

% lagged-spending instruments
Z = [X,Lagged,Lagged.^2];

% Gandhi-Houde nonlinear instruments
Z = GH_Instr([X(:,end),C,C.^2],Z,A);
Z0rfe = Z; % save for robustness checks
KZ = size(Z,2);

% first iteration
W = Z'*diag(xi.*xi)*Z; % GMM weighting matrix
Bblp = [Bols;0.25;0.15];
xi = xi_BLP(Bblp,xi,S,XC,A,nu,wnu);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit1_RFE] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% second iteration
W = Z' * diag(xi.*xi) * Z; % GMM weighting matrix
optVS=optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit2_RFE] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[B(1:length(Bblp),1);xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);
xi0rfe = xi;

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = -[XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z' * Ghat;
VBblp = inv(Ghat'*((Z'*diag(xi.*xi)*Z)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table 4 (column (V)):')
disp(print_results(round(rBblp_RFE,3),dnames,rfe,0,1,0))


%% region FE: BLP + neighbors' spending

% neighbors' spending and instrument
X = [X,[nb(:,1);nb(:,2);nb(1:N0,3);nb(N0+N1+(1:N2),3);nb(1:N0+N1,4);nb(:,5)]];
XC = [X,C,C.^2];
K = size(XC,2);
Z = [Z,[lnb(:,1);lnb(:,2);lnb(1:N0,3);lnb(N0+N1+(1:N2),3);lnb(1:N0+N1,4);lnb(:,5)]];
KZ = size(Z,2);

% BLP
W = Z' * diag(xi.*xi) * Z; % GMM weighting matrix
Bblp = [Bblp(1:K0);0;Bblp(K0+1:end)];
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_RFE_nb] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi;Z'*xi],[],[],...
    [sparse(KZ,length(Bblp)),Z',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = -[XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z'*Ghat;
VBblp = inv(Ghat'*((Z'*diag(xi.*xi)*Z)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE_nb = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table 4 (column (VI)):')
disp(print_results(round(rBblp_RFE_nb,3),dnames,rfe,1,1,0))


%%

save('voting_stage.mat')



